Direct Wigner tomography of a superconducting anharmonic oscillator 
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The analysis of wave-packet dynamics may be greatly simplified when viewed in phase-space. 
While harmonic oscillators are often used as a convenient platform to study wave-packets, arbitrary 
state preparation in these systems is more challenging. Here, we demonstrate a direct measurement 
of the Wigner distribution of complex photon states in an anharmonic oscillator - a superconducting 
phase circuit, biased in the small anharmonicity regime. We test our method on both non-classical 
states composed of two energy eigenstates and on the dynamics of a phase- locked wavepacket. This 
method requires a simple calibration, and is easily applicable in our system out to the fifth level. 



The state of open quantum systems is frequently rep- 
resented by a density matrix l2| • This description is 
useful for a wide variety of systems, in which the prob- 
abilities and coherences in some chosen basis are of in- 
terest. However, in systems with continuous degrees of 
freedom, e.g. the relative position and momentum of 
atoms in a molecule, wavepackets are often formed and 
therefore phase space distributions are better suited to 
characterize the state and its dynamics [3]. In particular, 
the Wigner distribution offers direct information about 
expectation values and purity, and provides a convenient 
framework to test the quantum-classical correspondence 
[4, 5]. In addition, since the Wigner distribution holds 
complete information about the state, it can be directly 
transformed into a density matrix and therefore its mea- 
surement is useful for quantum state tomography [6] as 
well. Numerous experiments measured the Wigner dis- 
tribution of harmonic systems using an auxiliary qubit 
[7, 8], or via heterodyne detection Q. While anhar- 
monic systems exhibit a wider variety of phenomena, a 
full quantum state reconstruction has so far been limited 
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In this letter, we report on a direct measurement of 
the Wigner distribution in the Josephson phase-circuit, 
a superconducting anharmonic oscillator. Our method 
utilizes simple tomography pulses, and as opposed to 
standard state tomography, does not require individual 
calibration of the pulses. 

Measuring the Wigner distribution in an anharmonic 
system poses several challenges. First, the phase of 
each level <p n in the rotating frame advances in an in- 
creasing rate with n, causing wavepackets to disperse 
in phase space during the tomography pulse. For ex- 
ample, in a cubic potential this rate is given by <p n « 
/3n(n — 1)/2, where /3 = 2tt (fio — /21) is the anharmonic- 
ity and n = 0, 1, 2, .... Second, in our system it is impos- 
sible to measure the probability distribution in the phase 
space quadratures directly, and therefore the prominent 
method for phase-space tomography is measurement of 
the parity after a coherent displacement [8]. In order to 
achieve an approximate displacement operation, one has 
to apply a pulse which is simultaneously resonant with all 



Figure 1: Experimental techniques, (a) The Josephson phase- 
circuit diagram, (b) Measurement of the occupation proba- 
bilities using short pulses in bias Imeas that cause selective 
tunneling of states n < k which are detected later with an 
on-chip SQUID, (c) Time sequence of the Wigner tomogra- 
phy measurement of Fock-type superpositions. 



the transitions within the measured subspace. We find 
that both restrictions can be practically met by both re- 
ducing the anharmonicity and applying sufficiently short 
tomography pulses. 

Our system, the Josephson phase circuit (see Fig. la), 
is a superconducting LC oscillator with a weak anhar- 
monicity, originating from a Josephson junction [14] . The 
potential energy of the system is one-dimensional and has 
the form of a double- well as a function of the phase differ- 
ence 5 across the junction. Using an external bias current 
lb we tune the anharmonicity and measure the occupa- 
tion probabilities of the energy levels inside the smaller 
potential well (see Fig. lb). At large anharmonicity 
(/?/27r/oi ~ 0.02), the lowest two levels are used for qubit 
operations via resonant microwave pulses I^ w that yield 
high-fidelity gates of only a few nanoseconds in duration 
[15j . At small anharmonicity (/3/2irfoi ~ 0.002), many 
more levels are accessible within a practical bandwidth 
and wave-packets can be easily generated [16] . In both 
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regimes the state is measured in the energy basis by ap- 
plying a short pulse in 1^ that causes selective tunneling 
of the probability in occupied levels n > k, where k de- 
pends on the height of the potential barrier, set by the 
pulse amplitude. Tunneling events are detected using 
an on-chip superconducting quantum interference device 
(SQUID), and for a set of measurement pulse amplitudes 
the occupation probabilities are reconstructed from the 
tunneling probabilities [16] . While standard state tomog- 
raphy (SST) can be used to measure the full density ma- 
trix of the lowest levels, in practice this method requires 
complicated calibration procedures for states with more 
than two levels, and is severely limited by decoherence. 
Our method requires only a single calibration measure- 
ment for all the tomography pulses, and can be scaled up 
to larger Hilbert spaces using optimized system parame- 
ters. 

Ideally, the Wigner distribution is proportional to the 
expectation value of the parity operator after a coherent 
displacement [17| . Due to the finite anharmonicity in our 
system, we use short, gaussian-shaped resonant pulses 
as an approximate coherent displacement, while work- 
ing at a weak anhramonicity. We control the phase and 
size of the complex displacement a by setting the phase 
and amplitude of the microwave pulse, using the rela- 
tion a = — (1/2) j Q(t)dt for a harmonic system, where 
Q(t) is the time-dependent, gaussian-shaped Rabi ampli- 
tude. For our experimental parameters (l8(] , the pulse is 
simultaneously resonant with many transitions between 
consecutive levels, ~ 6 of which are subject to an am- 
plitude variation of less than 10 % from the peak am- 
plitude [H]. For a given anharmonicity /?, a resonant 
pulse can be well approximated by a harmonic displace- 
ment in the limit j3T \ a\ m 2 /4 <C 1, where \a\ is the size 
of the applied displacement, m is the maximal occupied 
level after a displacement and T is the pulse duration 
(l9| . This condition limits the maximal displacement to 
be well below the size of the distribution in phase space 
(|c>W;e| <^ 1)? for initial states occupying up to 5 of the 
lowest levels. However, we find in simulation that while 
the phases of the displaced state are very sensitive to 
the above condition, the probabilities are not, and there- 
fore our approximation remains accurate for much larger 
displacements. We find that for our experiment parame- 
ters, up to 5 levels can be measured accurately. The ex- 
pectation value of the parity operator is calculated from 
the measured occupation probabilities and is given by, 

<* = ( 2 / 7r )E n (- 1 ) P n O- To test the effect Of this 

pulse on our system, we initiate our system in the ground 
state and measure the occupation probabilities immedi- 
ately after a short microwave pulse of total area a. Fig- 
ure 2 shows the results of this measurement as a func- 
tion of state number n and a (Fig. 2a), compared with 
the expected probabilities in a harmonic oscillator (Fig. 
2b), P(a, n) = (1/nl) exp(— \a\ 2 ) \a\ 2n . As expected, the 
probability distribution in n is narrower compared to the 
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Figure 2: Tomographic pulse. Occupation probabilities for 
levels n < 12, as a function of \a\ in the experiment (a) and 
theory (b). (c) and (d) show a histogram along the dashed 
lines shown in (a) and (b). 
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Table I: Results of the optimization algorithm. 



harmonic system for higher amplitudes due to our finite 
bandwidth. To compare the data and theory quantita- 
tively, we plot a histogram of the distribution (Fig. 2c, d) 
for a = —1.3 and a = 2.2. At a = —1.3 our data fits well 
with a harmonic displacement. At the largest displace- 
ment values (a = 2.2) the deviation from harmonicity 
becomes more apparent as expected [19] . 

To benchmark our method, we apply it on a set 
of superpositions of eigenstates of the system, = 
(|0) + e 1 ^ |/)) /y/2 [20]. These states have a simple struc- 
ture in phase space and cannot disperse since they con- 
tain one measurable phase that results in only a free rota- 
tion. As we operate our system in the small anharmonic- 
ity regime, preparing such states with high fidelity re- 
quires longer pulse sequences. Due to our short decay and 
coherence times (Xi = 120 ns, T 2 > 150 ns) we apply opti- 
mized control sequences to reduce decoherence effects to 
a minimum (see Fig. lc). We use a genetic optimization, 
where we define the population overlap y = Pideal • Pmeas 
and optimize it with a feed-back loop [19] , For a desired 
state |^), we set Pideal = (0.5, 0, .., 0.5, 0, ...) and use our 
measured occupation probabilities for Pmeas • The results 
of the optimization algorithm are summarized in Table 
1. 

As expected, x decreases for a larger n, mostly due 
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Figure 3: Superposition of Fock-type states, (a) Wigner tomography of genetically optimized superposition of fock-type states 
and (b) calculated Wigner distributions of these states, using level occupations from table 1 and phases of the extracted density 
matrices, (c) Extracted density matrices from the measurements shown in (a). 



to our relatively short preparation pulse duration r. To 
increase the value of x f° r these states we increase r, at 
the expense of the states' purity and optimization run- 
time. Our optimization algorithm does not take into ac- 
count phases, however for sufficiently high values of 
the Wigner image is dominated by only one phase, as 
can be seen in Fig. 3. In Fig 3a we plot the results of 
the tomography measurement on the ground state, and 
states described in Table 1. These are compared with 
the expected images (Fig. 3b), calculated using the mea- 
sured occupation probabilities, obtained independently 
by the optimization algorithm [21]. In Fig. 3c we plot 
the density matrices of the measured states, extracted 
using the harmonic oscillator wavefunctions. The diago- 
nal elements in the extracted density matrices agree well 
with the measured occupation probabilities up to level 
n = 4. However, simulation shows that off-diagonal ele- 
ments deviate increasingly for n > 2 [19]. We find that 
the origin of these errors is primarily the significant de- 
cay and decoherence that occur at highly excited lev- 
els during the tomography pulse, and secondarily the fi- 
nite bandwidth of our tomography pulse. The errors in 
the measured Wigner distribution can be further reduced 
with currently available higher coherence samples |19|. 

In general, for any finite set of repetitions, the amount 
of information that is gained about the density matrix is 
smaller in a Wigner measurement relative to SST (assum- 



ing equal number of repetitions and tomography pulses). 
This is due to the large overlap between the unitary op- 
erations associated with the Wigner tomography pulses, 
compared to the nearly orthogonal operations in a typical 
SST measurement. However, we find that the overhead 
in the number of repetitions in a Wigner measurement 
relative to SST with similar uncertainties is of order unity 
(l9| , thereby making the Wigner measurement a practical 
alternative for measuring the density matrix. 

Non-dispersive wavepackets are particularly interest- 
ing to measure in phase space, due to phase- locking be- 
tween the system and the drive that can take place within 
certain conditions [22] . In a system with negative an- 
harmonicity (/oi > /12)? a negative frequency-chirped 
drive (f drive < 0) can cause phase-locking at sufficiently 
large amplitudes, depending on the chirp rate and an- 
harmonicity. To measure this effect, we initiate the sys- 
tem in the ground state and apply a negative frequency 
chirp, with a drive amplitude above the phase- locking 
threshold and a final frequency centered close to the tran- 
sition frequency ^23- The chirp's temporal length and 
bandwidth \ ffi n — fin\ are chosen to be short (20 ns) and 
large (600 MHz) respectively, in order have a broad ex- 
citation of states and for the excitation to be adiabatic 
[~L6| . Figure 4 shows the result at selected times along 
the chirp (Fig. 4a), and after the drive is turned off 
(Fig. 4b). The axes in the images are rotated at each 



4 



10 ns 12 ns 16 ns 18 ns 20 ns 




-2-1 1 2-2-1 1 2-2-1 1 2-2-1 1 2-2-1 1 2 

Re(cc) 

10 ns 12 ns 16 ns 18 ns 20 ns 




Figure 4: Dynamics of a phase-locked wave-packet, (a) Wigner tomography during a chirp (T = 20 ns, Q/2ty = 66 MHz, 
fin — foi = 320 MHz, ffi n — foi = — 50 MHz), and (b) during drive-free evolution, (c) and (d) show the extracted density 
matrices for (a) and (b) respectively. 



time frame to fit the rotating frame of the drive. Dur- 
ing the chirp, we see a displacement of the ground state 
distribution that gradually acquires a constant phase as 
the drive crosses the linear resonance (/ = foi). This 
happens, as expected around t = 16 ns and the shape 
of the wavepacket becomes crescent-like. After the drive 
is turned off (t = 21ns), phase- locking is lost and the 
wavepacket disperses due to the finite anharmonicity. At 
t = 35 ns, the state has completely dispersed, however 
it still contains significant coherences, as indicated by 
the negative values in the Wigner plot and the large off- 
diagonal elements in the extracted density matrix (Fig. 
4d). The state then dephases into a ring shape which 
shrinks in radius towards the ground state. To track the 
decoherence dynamics in this experiment, we extract the 
state's purity evolution directly from the Wigner distri- 
bution: p = 7r J d~ot \W (~ot)\ 2 [3]. The result is shown 



in Fig. 5 (red circles), and compared with a simulation 
(solid line). As expected, the purity remains high dur- 
ing the chirp, and then decays as a result of decoherence. 
The purity reaches a minimum and then ascends towards 
unity with an exponential rate, consistent with the mea- 
sured decay time (Ti = 120 ns) and in accordance with 
simulation. 

In conclusion, we demonstrate a method of measur- 
ing the state of a multi-level system with a small anhar- 
monicity in phase space. Using this method we are able 
to accurately extract the density matrix of up to 5 levels, 
limited by systematic errors caused by decoherence and 
finite bandwidth. This represents a significant improve- 
ment on standard state tomography which is difficult to 
implement in this subspace. The tunability of the phase 
circuit offers an approach of improving the accuracy of 
such a measurement at larger anharmonicity, by quickly 
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Figure 5: State purity. Extracted state purity p from exper- 
iments (red circles) and simulation (solid line) during chirp 
and decay. The simulation includes the effect of the decay 
and dephasing and agrees with the measured decay time of 
Ti 120 ns. 

changing the bias 1^ after state preparation to the small 
anharmonicity regime where the tomography pulse is ap- 
proximately a harmonic displacement. Finally, we point 
out that the overhead in the number of measurements 
required to extract the density matrix in a Wigner mea- 
surement is relatively small, and in addition one gains 
global information about the state faster than in stan- 
dard state tomography. 

This work was supported by ISF grant 1248/10 and 
BSF grant 2008438. 
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Supplementary Information 

Materials and methods. The Josephson phase cir- 
cuit [lfl used in the experiment has the following design 
parameters: critical current Io « 1.5 /iA, capacitance 
C~1.3pF and inductance L^940pH. The qubit has a 
tunable frequency /oi in the 6-9 GHz range. During the 
experiment the device is thermally anchored to the mix- 
ing chamber of a dilution refrigerator at 30 mK, where 
thermal excitations of the qubit are negligible. 

We use a custom built arbitrary waveform genera- 
tor (AWG) having a fast (Ins time resolution), 14-bit 
digital-to-analog converter to shape the microwave pulses 
which control the quantum state of the circuit. We con- 
trol the phase and amplitude of the drive by modulat- 
ing a high-frequency oscillator of frequency /lo with 
the AWG, using an IQ-mixer. The modulation sig- 
nals are fed into the I and Q ports of the IQ-mixer to 
give an amplitude ^I(t) 2 + Q(t) 2 and a relative phase 
4>(t) = arctan (Q(t)/I(t)) at its output. 

Harmonic response. As indicated in the 
manuscript, the bandwidth of the tomography pulse must 
be large relative to the anharmonicity, for a nearly har- 
monic response. This is easily understood by looking at 
the position of the transitions in the spectral domain, 
relative to the local oscillator frequency (see Fig. 1ST]) . 
As seen in the figure, the effective drive amplitude of the 
lowest 6 transitions is the same, to within a 10 % vari- 
ation. This is consistent with our data and simulations 
(see Sec. II), showing small systematic errors within the 
corresponding subspace. In principle, more transitions 
can fit within this region by positioning the local oscilla- 
tor lower in frequency. 
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Figure SI: Tomography pulse envelope in the frequency do- 
main. Dashed lines are the lowest five transitions f n ,n+i and 
the local oscillator frequency position in the experiment. The 
pulse envelope (solid line) is the normalized fourier transform 
of a 1.6 ns FWHM gaussian function. 




Figure S2: Measurement pulse used in the experiment. All 
the pulse parameters except for the total amplitude A are 
fixed. 



Populations measurement. The measured escape 
probabilities due to a set of experiments with different 
measurement pulses are converted into level occupation 
probabilities using the measured escape curves of the low- 
est levels in the well [2|. In this experiment, we used a 
novel measurement pulse shape, to reduce distortions of 
the escape curves of higher levels. This measurement 
pulse is shown in Fig. [S2l Typically, one must use a 
pulse with a slowly decaying end, in order to reduce the 
effect of population re-trapping in the original well af- 
ter the pulse ends. The effect of re-trapping becomes 
more prominent as the anharmonicity is reduced due to 
the corresponding deepening of the potential well, and 
therefore requires increasing both T2 and T3 to eliminate 
re-trapping. At the working bias point used in our exper- 
iment, we use T2 = 25 ns and T3 = 15 ns. These param- 
eters are non-negligible relative to the decay time and 
cause artificial increase of the extracted population at 
lower levels, which in turn distorts the Wigner image. 
To avoid this effect we add a small step in amplitude at 
the beginning of the pulse (r± = 1 ns) to preselect the 
population we want to escape, and then reduce the am- 
plitude by an amount 8 A. Even the small difference 8 A is 
sufficient to protect undesired population at lower levels 
from tunneling out during the long waiting time in which 
the escaped population decays in the external well. 

Density matrix fit. In the Wigner tomography ex- 
periments we extract the density matrix from the pop- 
ulations of the displaced states [3|. We use 200 ho- 
mogeneously distributed random displacements within a 
\a\ < 2 circle to fit the density matrix, while restrict- 
ing the density matrix to a 6 x 6 subspace. It should be 
noted that both the measured Wigner distribution and 
the extracted density matrix represent the state after a 
rotation that occurs during the tomography pulse. To get 
more accurate phases, one can apply an inverse propa- 
gator on the density matrix U = exp(—iHoAt/h) : where 
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Hq is the drive-free Hamiltonian and At is the effective 
pulse length for the rotation. 



STATE PREPARATION BY GENETIC 
OPTIMIZATION 

Finite anharmonicity makes it possible to prepare 
states, composed of arbitrary superpositions of eigen- 
states within our system. However, fourier broadening 
of the drive, together with power broadening of the tran- 
sition energies causes nontrivial excitation at relatively 
small anharmonicities. Adding the short decay and co- 
herence time, it becomes challenging to prepare a de- 
sired, yet pure state. We solve this difficulty, by op- 
timizing the state produced with a feedback from the 
experiment. Our target state in the optimization algo- 
rithm is a superposition of eigenstates of our system: 
= (|0> + e i(t> |0) /y/2. Using the measured proba- 
bilities we are able to optimize such a superposition up 
to an unknown phase <p that causes only a free rotation 
of the measured Wigner function. 

Our optimization algorithm is based on guided evolu- 
tion. It is constructed in the following steps: 

(a) We define a set of Nq pulse sequences (genomes). 
Each genome is a sequence of N t complex numbers that 
represent the amplitude and phase of the drive's pulse 
envelope at each time step, with 1 ns time resolution. In 

addition, each genome is associated with a population 

> > — - — > 

overla p, defined as: x = Pideal ' Pmeas, where P ide ai and 

Pmeas are the desired and measured population vectors . 
The genomes are initialized with random complex num- 
bers, each having a maximal amplitude £l m axi and x — 0. 

(b) At each iteration of the algorithm (generation) , we 
assume that the current set of genomes is sorted in order 
of decreasing x- The upper N a genomes are kept, and 
the next N a genomes are replaced by a combination of 
the upper N a + 1 genomes: G Na +l = ~^(G k , Gfc+i) + 
where the function ~P(x,y) randomly selects an ampli- 
tude between the sets x and y at each time step and ~^ 
is a set of random complex numbers with a small ampli- 
tude e max (noise). The bottom Nq — 27V a genomes are 
replaced by random complex numbers with an amplitude 

^max • 

(c) The bottom Nq — N a genomes of the next gen- 
eration are then applied and the resulting probabilities 

Pmeas f° r eacn genome are measured. After calculating 
X of each, we sort our new set G in order of decreasing 
X, and mark genome sequences that have (potentially) 
higher x than the current optimal x- For each genome 
having a potentially higher x we repeat the measurement 
N rep times, recalculate x and reposition them in the set 
G. Step (b) is run again. 

Fidelity measure. Our population overlap x mea- 
sures the distance between two states. There are many 



possibilities of defining the fidelity, and in general each 
gives a different bias for the algorithm towards finding a 
specific family of states. In our definition for example, 
X is highly sensitive to the overall population that ex- 
ists outside the subspace that is populated in the target 
state. However it is much less sensitive to the distribu- 
tion among levels that are populated in the target state. 
Our definition is particularly useful for optimizing our 
target states \iJji) for a Wigner measurement, because we 
wish to eliminate the population of states other than |0) 
and \l). 

There are two important factors that limit the per- 
formance of any optimization algorithm in experimental 
systems: shot noise, and drifts in the physical parameters 
of the system. 

Shot noise. Any finite number of repetitions results 
in some uncertainty in the measured probabilities. To de- 
termine n occupation probabilities in the final state we 
measure the tunneling probabilities after n measurement 
pulses having different amplitudes. In a typical optimiza- 
tion algorithm, we repeat the measurement of a single 
tunneling probability 900 times. Statistical analysis gives 
a typical uncertainty (standard deviation) of ~2 % in x 
for sufficiently high values of x (x ^ 80%) and there- 
fore one cannot distinguish between two genomes with a 
X difference which is smaller than ~3 %. We therefore 
expect the algorithm to run significantly more slowly at 
X's approaching unity, due to false increase of x- To in- 
crease the efficiency of the algorithm without increasing 
the number of repetitions of each measurement, we re- 
peat the measurement only for genomes with potentially 
increased x? as described in step (c). 

Drifts. The response of the system to the application 
of a given genome can change in time due to drifts in its 
physical parameters. We correct for drifts in the energies 
by performing a spectroscopic measurement of the qubit 
transition frequency every 10 minutes and adjusting the 
flux bias to restore the qubit frequency to its original 
value. In addition, we correct for drifts in the offset volt- 
ages applied to the IQ mixer that eliminate leakage of 
the microwave drive while it is turned off. 

SYSTEMATIC ERRORS 

To quantify the errors in the measured Wigner function 
and density matrix caused by finite anharmonicity and 
decoherence, we perform numerical simulations . In our 
simulations, we propagate an anharmonic system, initial- 
ized with a pure state p = \ip) (ip\, with resonant pulses. 
Each pulse is assumed to be resonant with the first tran- 
sition and has a gaussian envelope with the same width 
as in the experiment. We use the rotating wave approx- 
imation for constructing the Hamiltonian, and neglect 
corrections to the drive coupling beyond the harmonic 
approximation, namely (n\S\m) = for n ^ m ± 1 and 
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Figure S3: Systematic errors in the Wigner distribution, for 
initial states \ipi) = (|0) + |/)) /y/2, where I = 1,2,3,4. The 
cross-correlation deviation as a function of the anharmonicity 
(a) and initial state. In (a), decoherence is not included, and 
the magenta, red, blue, green lines correspond to the initial 
state I = 1, 2, 3, 4 respectively. In (b) decoherence is included 
and the anharmonicity is fixed to /3 — 20 MHz. The black and 
grey lines correspond to decoherence parameters Ti — 150 ns, 
T 2 = 200 ns and Ti = 600 ns, T 2 = 600 ns respectively. 



(n\ 5 \m) = ^Jri, y/n + 1 for m = nil. This results in 
the following time dependent unitary propagator: XJ(t) = 
exp (if [(Q(t)a^ + Sl(t)*a) + fia^a(a^a - 1)] ) , where St 
is the time step in the simulation, Q(t) is the time de- 
pendent drive envelope amplitude and /3 = 2n (/oi — /12) 
is the anharmonicity. Decoherence and energy relaxation 
are taken into account using quantum operations, assum- 
ing only two phenomenological parameters: the energy 
relaxation time of the qubit Ti and its pure dephasing 
time T 2 . 

We expect to have negligible systematic errors due to 
finite anharmonicity when the second term in the ex- 
ponent of the propagator becomes negligible. By sep- 
arating the terms in the exponent to first order, and 
assuming a constant drive amplitude Q of total dura- 
tion T, we can approximate the propagator to U ~ 
D(a) exp(— i^-UaaJ — a*a) , a^a(a^a — 1)]) where D(a) 
is the displacement operator and a = iQT. The sec- 
ond term can be considered negligible in the limit, 
\a\Tf3m 2 /A <C 1, where m is the maximal occupied 
state. For our experimental parameters (P/2tt = 20 MHz, 
T = 1.6 ns), the second term can be neglected only for 
\a\ <C 1, however, the error in the Wigner distribution, 
obtained from the state populations after the pulse are 
negligible even for a « 2, as we show in our analysis. 

We plot the errors in the Wigner distribution and the 
extracted density matrices, as a function of the anhar- 
monicity /?, initial state and decoherence parameters. 
The density matrix is extracted, as in the experimental 
data, using only the populations of the lowest 6 levels. 

To quantify the error in the Wigner distribution, 
we calculate the cross-correlation (at zero offsets) be- 
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Figure S4: Systematic errors in the extracted density matrix, 
for initial states \ipi) = (|0) + |/)) /y/2, where I = 1,2,3,4. 
The left column of plots shows the fidelity error, and errors 
in the off-diagonal matrix element poi as a function of anhar- 
monicity, in the absence of decoherence. The magenta, red, 
blue and green lines correspond to the initial state / = 1, 2, 3, 4 
respectively. The same error measures are plotted in the right 
column, as a function of initial state /, with decoherence in- 
cluded, at P = 20 MHz. The black and grey lines correspond 
to decoherence parameters T\ — 150 ns, T2 = 200 ns and 
Ti = 600 ns, T 2 = 600 ns respectively. 



tween the ideal Wigner distribution of the initial state, 
and the one obtained from the expectation value of 
the parity operator after a set of displacements. The 
cross correlation is defined as C(f(x,y),g(x,y)) — 

(f(x',y')-(f))(g(x' jy ')-(g)) 



E 



E (/(*" ',y")-a» 2 E(p(*' y \y")-{g)) 2 



where (f),(g) 



are the average values of /, g. The results are plotted in 
Fig. [SU As seen in the figure, the error increases sharply 
with both anharmonicity and maximal populated level I. 

To quantify the error in the extracted density ma- 
trices, we use the standard fidelity definition: (SJ 7 ) = 

1 — Tr\J \/ r p^ p 1 \/~p^ , where p 1 is the initial density ma- 
trix and p F is the density matrix obtained from a fit to 
the populations of the displaced states. We choose the 
initial density matrix p 1 to be |^)(^|, where = 
(|0) + \l))/y/2. In addition, we calculate the error in 
the non-zero off-diagonal elements of the density matrix 
using the following definitions: Spoi = I Ipozl ~ |/°oz| I ^ s 
the error in the amplitude of the matrix element, and 
S(f)(poi) = (pQi) — <fr (poi) I i s ^ ne normalized error 
in the phase of the matrix element. 

As seen in the Fig. |S4} all the error measures are 
negligible (smaller than ~ 0.05) when decoherence is 
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Figure S5: Error in the fidelity due to shot noise in standard 
state tomography and WT. (6 J 7 ) vs. N in WT (blue) and 
SST (red) of the initial states = (|0) + |4)) /V2 (a) and 

|^) = (I/a/5) E l&) ( b )- Both follow the VV^V trend. 

neglected. However, when included, both the fidelity 
and the amplitude of the off-diagonal elements have non- 
negligible errors. For currently available samples having 
T\ > 600 ns, and correspondingly longer T2, the errors 
due to decoherence can be substantially reduced; In this 
case, the errors become smaller than 0.1 in all the mea- 
sures. As expected, the phase error is quite insensitive 
to decoherence. 

SHOT NOISE IN WIGNER TOMOGRAPHY VS. 
STANDARD STATE TOMOGRAPHY 

As pointed out in the manuscript, phase space is a 
convenient basis to experimentally acquire global infor- 
mation about the state (e.g. phase distribution, average 
energy etc.) fast, on the expense of accurate knowledge 
of the state in the eigenstate basis. Therefore, when ex- 
tracting the density matrix in the latter basis from a 
Wigner measurement, one requires an excess number of 
measurements compared to standard state tomography 
(SST) to achieve similar uncertainties. It turns out that 
this is only true for a certain class of states. For states 
that are dispersed in phase space, such as states com- 
posed of a coherent superposition of a small number of 
eigenstates, SST requires significantly less tomography 
pulses than Wigner tomography (WT) to achieve a com- 
parable error in the density matrix but the same num- 
ber of pulses for states that are localized in phase space. 
The overhead in the number of measurements is an im- 
portant parameter from an experimental standpoint, and 
therefore we perform numerical simulations to calculate 
it. In the following we describe our simulation methods 
for both cases. 

WT. We start with a pure initial state We use 
the same method described in Sec. II to calculate the 
density matrix after a coherent displacement. We keep 
the diagonal elements of the final density matrix, for a set 
of Nw random displacements a that are uniformly dis- 



tributed in the complex plane. For each experiment (a 
particular displacement) we calculate an ensemble of M 
possible outcomes for the measurements of the diagonal 
elements, assuming r repetitions in the experiment, and 
a binomial distribution for the escape probabilities, from 
which the diagonal elements P(n) are calculated. For 
each outcome, we extract the density matrix and calcu- 
late its fidelity T = Tr^J v /r p F ^ p ideal y^p^ . We then 
calculate the average fidelity (J 7 ) of the ensemble to find 
the experimental error (8 J 7 ) = 1 — (J 7 ) due to shot noise. 

SST. We start with a pure initial state We con- 
struct a set of Nsst = d 2 orthogonal, unitary operations 
Uj to span a d-level subspace. The set of operations is 
chosen to be the generators of SU(n) for convenience. 
From the diagonal elements of the rotated density matri- 
ces, we are able to extract the expectation values of each 
generator (Uj) = Tr(pUj), and therefore reconstruct the 
original density matrix: p = ^2(Uj)XJj [4|]. As before, 

3 

for each diagonal we calculate an ensemble of M possi- 
ble measurement outcomes due to shot noise. From each 
of the measurements in the ensemble, we extract the ex- 
pectation values of the operators Uj and calculate the 
corresponding density matrix. The average fidelity error 
of the ensemble of density matrices is then evaluated. 

Figure [S5] shows the results of both simulations. We 
plot (5 J 7 ) as a function of N = R, where R is the number 
of repetitions of the experiment per tomography pulse. In 
the WT simulation we define N relative to the SST case: 
N = R (^ n^ t ^ ? where Nw is the number of tomography 
pulses in the Wigner simulation, and Nsst is the number 
of tomography pulses in SST. While Nsst is fixed, we 
vary N in the WT simulation by using many displace- 
ment pulses while fixing R, and in the SST simulation 
we vary TV by changing the statistics R. In all Wigner 
simulations we fix R to be 900 (as in our experiment) and 
in the SST simulations we vary R from 100 to 3000. 

We see that for an initial state, composed of a coherent 
superposition of only 2 states (Fig. IS5k), SST outper- 
forms WT by a factor of 8, in terms of the fidelity J 7 . 
In contrast, for states composed of a coherent superposi- 
tion of 5 eigenstates, the amount of information that is 
extracted per pulse in SST and WT is similar. This is 
because our chosen state is partially localized in phase 
space, and therefore requires less displacement pulses to 
extract its properties. 
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